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Abstract 

We introduce a new formalism for dealing with networks of queues. The formalism is based on the 
Doi-Peliti second quantization method for reaction diffusion systems. As a demonstration of the method's 
utility we compute perturbatively the different time busy-busy correlations between two servers in a Jackson 
network. 
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1 Introduction 

A Jackson network consists of a set of customers performing random walks within a network of queues, where 
customers may also be added (in a Poisson manner) or removed from the network. A well known theorem 
(Jackson's theorempQ) states that the instantaneous steady-state behavior of such a queueing network has 
the same statistics as a set of independent M/M/l queues. In contrast to this simple result, dynamical 
properties of queueing networks are quite involved due to the appearance of correlations between events at 
different times. As very little is known about the nature of these correlations, a systematic approach would 
be interesting. Providing such an approach is one of the purposes of this paper. 

Another purpose of this paper is to point out a connection between queueing networks and various physical 
models known as "reaction diffusion models" (RD) that are used to model bulk chemical reactions and various 
other particle dynamics. RD systems have benefited enormously from a major insight into the problem that 
was made by Doi and Peliti, who noticed ( independently) the usefulness of quantum many body techniques 
to analyze the RD problem. Their insight was that the RD problem could be written using quantum- 
mechanics-like "second quantized" operators to describe the hopping and interactions of the particles. 

RD models are used to describe the microscopic motion of particles through a medium which has a diffusive 
effect on the particles. If one considers a queue customer as a particle and a server as a site at which an 
interaction takes place, the hopping can be considered as the result of the customers getting randomly 
routed to other servers and being queued there. The description of a queueing network as a set of customers 
performing random walks is more or less standard, but the relationship to molecules adrift in a medium 
seems to have been neglected. Once the analogy is made, we can borrow some of the techniques used in 
these theories to reformulate some queueing models in a suggestive (and in some cases simpler) form. This 
reformulation will allow us to develop a systematic perturbation expansion of server state correlations i.e. 
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the busy-busy correlation of two different queues (at different times) . We shall calculate the first order term 
in the perturbation, which appears to be a new result, and gives us a formula that is found to be numerically 
valid across a wide range of simulated network parameters. 

The paper is organized as follows: after a brief review of relevant previous work the Doi-Peliti formalism 
shall be introduced and applied to queues in section [3l In sub-section 13.31 we shall introduce the dynamic 
operator that represents the Jackson network and show how the equilibrium statistics can be easily derived 
from the operator. The generalization to dynamic correlations will be discussed in section |H In section [5] we 
shall generate a formal perturbation expansion for the "propagator" of the Jackson network operator. We 
shall use the expansion to calculate the (Laplace transformed) "busy-busy" correlation function. In section 
[6] the perturbation results will be compared to simulations. Section [7| is devoted to remarks on possible 
extensions and shortcomings of the technique. 



2 Previous Work 

Second quantization was introduced in reaction diffusion problems by Doi and Peliti|2j. The method has 
been developed quite extensively by Cardy and othersp]. The relationship between RD and queueing was 
pointed out by [4]. 

The possibility to compute correlation functions perturbatively for the queue network depends on the 
knowledge of the Green's function for the single M/M/l queue. A convenient representation, that maps into 
a normal-ordered second quantized form was derived by [7]. 

An approach that is very similar to the one presented here is that of Massey[5j[6j. In a series of papers 
he defines "an operator theoretic approach" to Markovian queues. These operators are reminiscent of the 
quantum mechanical (QM) creation and annihilation operators introduced here. Indeed, the work presented 
here could be viewed as an extension of Massey's work (although the author came about the representation 
independently), but the emphasis of this paper is quite different. This paper deals mostly with dynamical 
properties of queues and the possibility to get perturbative results for correlations. While section [3] can be 
viewed as a review of known queueing results in a second quantized framework, the rest of the paper goes 
beyond this, to obtain new, dynamical results. For readers unfamiliar with the operator theoretic approach, 
section [3] can be used to make contact with standard queueing theory formalism. 



3 Second quantization formalism 
3.1 General 

We define the state of a queue as the number of particles (customers) stored in the queue (we include the 
customer being serviced as being the first in queue). A queue with n particles is denoted by the "bet" \n). 
In addition, we define a set of orthogonal "bra" states (m| such that the inner product (m|n) = S mn • For a 
system of multiple queues indexed by i a snapshot of the system is given by a direct product of all single 
site states: |m, n-z . . .) — Y\ k \n k ) k . The probability vector of all states can be given as a sum of all possible 
configurations \ip) = J2 ni Em • ■ • p ( n u n 2 • • • ) life K)- 

Following the usual course[33, we define creation and annihilation operators (sometimes called " ladder" 
operators) a + and a respectively, that have the following effect on the states: 
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a\n) = n\n — 1) annihilation 
a + \n) = \n + 1) creation 



(3.1) 



Which generalizes to 



ai\n 1 ) 1 \n 2 }2 ■ ■ ■ |n<)< ■•• = n;|ni)i|n 2 )2 ... |n< - 1), . . . annihilation 
a+|m>i|ri2>2 • ■ ■ K)i ' ' ' = |ni)i|n 2 ) 2 • ■ ■ + l)i . . . creation (3.2) 

It turns out that although the annihilation operator a has a simple commutation relation with the creation 
operator a + , it is not a useful operator when dealing with the M/M/l queue. This is due to the fact that he 
probability to leave a queue is not a function of the number of items in the queue (no mass action law[8j). 
Instead of the a operator, we shall define an modified annihilation operator Q such that 



[n — 1) = Q\n) for n > 

= Q|0) (3.3) 

So that a + can be considered as an operator that adds a single client to a queue and Q as an operator 
that removes a client (i.e. via serving the client). 

3.2 Single queue M/M/l at equilibrium 

Using the ladder operators Q and a + [9] , the M/M/l queue master equation 



P k = n(P k+ i - P k ) + A(P fc -i - P k ) k>0 

P = fxPx - \P (3.4) 



has the operator form of 



dt\ip)=£\ip) (3.5) 

with 

C = (1 - a+)(nQ - A) = /i(l - a + )(Q - p). (3.6) 

Where p — as usual. 

There is nothing "mystical" about the notation. The summation of configurations that define 
the "wave function" is quite similar to a z-transform of the probability vector, i.e. ip(z) = 
Yl n Sn, ■ ■ ■ -P( n i: n 2 ■ ■ • ) Ilfc z k k ■ The ladder operator a + is equivalent to multiplying by z and the lowering 
operator a is equivalent to differentiating by z (that is d z ). The Q operator is slightly less familiar, but 
can be viewed as the operation of ip{z) — > 3^MzJ^2] _ jg j us t more convenient to treat these as abstract 
operators. 
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As pointed out in Massey's papers, the equilibrium solution 



= C\if>) (3.7) 

can be obtained by observing that in order to obtain = /i(l — a + )(Q — p)\ip) it is enough to find a state 
that obeys = (Q — p)\ip) ■ We note that the modified annihilation operator Q generates a set of eigen-states, 
similar to the "coherent states" that exist in quantum mechanics. Namely 

|0) =Q_L_|0>. (3.8) 



1 - xa+ 1 ' 1 

As can be seen by expanding in x. We see that (Q — p) l _ 1 xa + |0) = (x — p) 1 _ 1 xa + |0) so that by setting 
x = p we can solve equation 13.71 with 

\P) = T^^IO) = {1 - P}{|0) + p\l) + P 2 \2) + . . . } (3.9) 
1 — pa + 

The formal time dependent solution to the M/M/l queue, starting with an initial probability state \tjjo) 
is \ip(t)) — e tc \ipo) and due to the Markov-chain nature of the M/M/l model we expect that the long time 
behavior of arbitrary physical initial conditions relaxes to the equilibrium state |/o) [lO] - In what follows we 
shall sometimes call the time domain Green's function e tc the propagator. Various average quantities can 
be obtained by constructing expectation values with the "unit bra" 

(J| = (0| + <1| + (2|+... (3.10) 

For example, suppose we are interested in the average busy ratio b of an equilibrized queue. The operator 
combination a + Q represents an object that returns 1 when the queue is not empty and otherwise. Hence 
b = (I\a + Q\p) which can be easily evaluated as p . 

3.3 Open Jackson network of queues at equilibrium 

The open Jackson network looks like this: 

d t W = 5^(1 - af)[J2(Sij ~ U^iQi - 7j M (3.11) 

3 i 

To simplify, we define Lij = (Sij — r^j)pi and 

LijPt = Jj (3.12) 
so that the open Jackson network operator becomes 

dt\il>) = (1 - a^LijiQi - piM (3.13) 
In this form is now quite easy to verify that the product form 

\pi, Pi-, ■ ■ ■ ) = life 1^+^ |0) is the stationary solution of the open M/M/l Jackson network, thus proving 

1 a k P k 

Jackson's theorem. 



4 



4 Queue dynamics in the second quantized formalism 



4.1 Expectation values and correlations 

We have seen that the second quantized formalism allows us to express averages of functions of the occupation 
number as the expectation values of various operators. For example, as stated above, the operator a + Q 
projects out of a probability vector all states that are non empty (namely a + Q\n) = |n) for n > and 
a + Q\0) = 0). Similarly, (a + ) 2 (Q) 2 projects out all states that have 2 or more customers in queue. In 
fact, the term (a + ) n Q n — (a + )" +1 Q n+1 projects out the state that has exactly n customers. We see that at 
least formally, given a probability vector \ip)we can extract all queue occupation information by considering 
expectation values of operators O, (l\0\ip). 

Suppose we know the queue state at time and we want to find the average of some queue occupation 
number quantity at time t . Assuming that the initial queue state is |^o)j the formal solution at time t is 
\ip(t)) — e 1 \ipo) so that the expectation value of the operator O becomes (I\0\ip(t)) = (I\Oe tc \ipo) . This 
expectation value can be interpreted as an average measurement performed on a queue after it has evolved 
from the initial state for a period t. 

Let us now consider the case of correlation functions, i.e. measurements taken at two different times. If 
we have two operators 0\,0-2 that represent two measurements of the queue behavior, we can construct the 
correlation function as (I\02e tL Oi\ip) . This object is equivalent to starting out with a distribution of 
initial conditions of the queue, measuring the value of 0%, then evolving the queue for a time t, measuring 
the value of O2 and then averaging over all evolutions and all initial conditions given by 1^). If we are 
interested in steady state correlations, we may replace the generic initial condition with the steady state 
solution \p) . 

Correlation functions are important because they give dynamical information regarding queue evolution. 
The correlation function is the first moment of the joint probability distribution < ab >= J dadb abP(a, b), 
so that the joint probability of our two measurements can be extracted from a generating function J(p, q) = 
(I\e lp02 e tL e lq02 \ip). 



4.2 Single queue M/M/l correlation functions 

What is the correlation between the server state (i.e. busy or empty) at time and time t ? If the initial 
queue state distribution is |-0o) we need to compute 

_ (I\a+Qe tc a + Q\^ ) 

C[t) ~ (ibM 

which simplifies for |^o) = \p) to C{t) = p(I\Qe tc a + \p) . 

Quantum mechanics teaches us that expressions such as this can be dealt with conveniently if the operators 
inside are "normal-ordered", that is, brought to a form such that all annihilation operators are near the ket 
\p) and all the creation operators are near the bra (I\ . The propagator g(t) = e tc does not have a simple 
normal ordered form, but happily, its Laplace transform g(uj) = -^zrr, nas been reduced to such a form by [7], 
although they did not use an operator formalism. For completeness we present the derivation in appendix 
lAl The final result is 

gfcj) = X /^P\ (i + - a+Q))- 1 —— (4.2) 

yy ' (l-xa+) y p-x K ^"{1-xQ/p) y ' 

where 
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x(uS) 



(u/n + p + 1) - y/jpfji + p + l) 2 - 4p 
2 



(4.3) 



We note that g(uj) is normal ordered and that computing C(ui) 
(although tedious) exercise if one recalls that Qa + — 1. The result is 



p(I\Qg(oj)a + \p) becomes a simple 



p{I\Qg(u)a+\p) 



x 



(p + px 



x) 



(4.4) 



M(l -x)(p-x) 



(see appendix [B] for details). 
4.3 Jackson network correlation functions 

Consider queue (3 at time and queue a at time t (a^ If both queues are not empty at the respective 
times, we say the busy-busy correlation is 1, and zero otherwise. In equilibrium ( that is, starting out with 
the stationary state at time zero), this function can be described as: 



the exact normal ordered form for e or for its Laplace transform ^jzrjr is n °t known, so that in order to 
calculate the correlation function we must resort to approximations. In the next section we shall apply 
perturbation theory, which allows a systematic approximation scheme in terms of a small parameter. 

5 Perturbation expansion for Jackson networks 
5.1 The General Formalism 

We shall consider a perturbative expansion around the diagonal part of the Jackson operator. That is, we 
consider the Jackson network as a perturbation around a set of independent M/M/l queues, with a weak 
coupling to each other. The source rate of each unperturbed queue is pkPk, that is, the effective rate which 
appears in the steady state solution to the full network. The perturbation is the non Markovian effect of 
cross talk between the independent queues, 
if we consider 



where this time, C = (1 — a 




Lij =L° + eL 1 = (Sij)pi 



er, 



(5.1) 



we may write the operator as 





(5.2) 



L" + eW 



6 



Where we have set L° = 52j(l ~ a ^)^j[Qj ~ Pj] as the unperturbed operator and W = — Ylijfi ~ 
a1j )ri^j^i[Qi — pj] the perturbation. 

While it may seem unnatural to fix the values of the pk and to expand around the diagonal terms of L 
rather than fixing the values of the 7fc, this form of splitting of the operator has the advantage that the 
stationary solution remains unchanged for all values of e. When we calculate correlations with respect to 
the stationary state, we can avoid modifications (familiar in perturbation theoretic expansions) due to the 
change in the stationary state. 

The propagator for the unperturbed operator is just a product of propagators for uncoupled queues 

0°(*)=nsi°(t) 

The propagator e tc becomes 

eW+'W =g°(t)+e [ dTg°{T)wg°(t-T) + e 2 [ d n r dT 2 g°(T 1 )wg {T 2 )wg°{t-T 1 -T 2 ) + ... (5.3) 

Jo Jo Jo 

and Laplace transforming yields: 

G(u) = G°(w) + eG°(w)WG°(w) + e 2 G°{cj)WG°(uj)WG {u) + ... (5.4) 

While trying to obtain useful perturbative results we are faced with two technical difficulties, the first 
being the need to put the terms in the perturbation series in normal order and the second is the need to 
Laplace transform a product of propagators. The transformation of the product form of Q {t) — Yidi (*) 
into the convolved Laplace transformed form of G°(ui) is rather involved. Dealing only with off-diagonal 
correlations and only the first order in the perturbation expansion will allow us to simplify matters, as long 
as we consider direct calculations of correlation functions. 

5.2 Perturbation expansion of the busy-busy correlation function 

Consider queue (3 at time and queue a at time t (a ^ /?). If both queues are not empty at the respective 
times, we say the busy-busy correlation is 1, and zero otherwise. As pointed out in 14.31 the busy busy 
correlation function is: 

C a0 (t) = (I\a+Q a g(t)a+Q \p) = p (I a \Q a g(t)a+\p a ) (5.5) 

Expanding C a ^{t)'m e as C a ^(t) — {t) + eC^ {t) + . . . we find (using the propagator expansion in 
the time domain (equation (|5.3|l ). 

Cl p {t) = p (I a \Q a g°(t)a+\ Pa ) 

C^(t) = p [ (I\Q a g Q (T)Wg°(t-T)a+\p) (5.6) 
Jo 

Furthermore, the product structures of g° and of the stationary state \p) allow us to decouple the zero 
order term in the expansion: 

C%(t) = pp(I\Q a g o a (t)g° (t)a+\p) = p p (I\Q a a+\p) = Pp <Q a >< a+ > (5.7) 
and the first order perturbation is: 
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Ct^wE^J I (IlQag^il-apiQi-p^it-^a+^dr. (5.8 
All cases where a ^ jor j3 ^ i vanish due to the commutations [X,-, Yj] = for i ^ j, leaving us with 



C* /?(*) = PpSL p , a / (I\Q a g° a (T)(l - a+)[Q - p p ]g° (t - r)a+\p)dr. (5.9) 
Jo 

Again, the product structure of the stationary solution comes to our aid allowing us to decouple: 

Cl p {t)=p 5Lp, a f (I a \Q a gl{T)(l-ai)\p a )(Ip\[Qp~pp]gl(t^T)a + \pp)dr. 
Jo 

This is simple enough to allow us to perform a Laplace transform: 

C af) {u)=C%{u) + eC 1 atP {w) + .... 

with 



AO / \ _ Pi3(I\Q a ap\p) _ 

ot,p\ / UJ .. 

= P/3^/3^< / Q |Qc e 5a(w)(l-a+)|p a )(^|[g^-p^]^(t i ;)a^|p / 3) 
Expanding the C7* Juj) expression we remain with 

C afi (w) = + epp6Lp, a [— - {I a \Q a g°(u))a+\p a )][(I \Q g° (u))a^\p ) 



^1+eppSLpA- ~ (I a \Qa^)a+\p a )][(Ip\Qpf p (u)a+\pp) - ^] 

UJ UJ UJ 

and using the results in appendix iBl 



C a p(u)) = ^ + epp8L^ a [ -}[ - + ] + ... (5.10 

UJ UJ PaPa{l-X a ) ppppil ~ X{3) U) 

This is the main result of the paper. The non zero O(e) term in the expansion is a direct demonstration 
of the fact that the actual dynamics of Jackson networks are not equivalent to an independent Poissonian 
arrival system, even though the stationary state appears to behave as one. In the next section we shall 
compare our result to simulations of various networks. But before that we shall comment on the correlation 
behavior a single queue with respect to itself. 

5.3 Same queue correlations and busy periods 

The same queue correlation function 

C aa (t) = (I\aiQ a G{t)a+Q a \p) = p a (I\QaG(t)ai \p) cannot be deduced as special case of the inter-queue 
correlation calculated above. This is due to two facts. The first is that the unperturbed situation C® a (t) is 
different: C® a (t) — p a {I\Q a g^ (t)a+||o) , depending on one g° only. The second fact is that the first order 
term in e vanishes because the perturbation term W doesn't contribute to the diagonal, W aa = 0. 

It is beyond the scope of this paper to deal with the second order contribution, however, an interesting 
result can be extracted from this behavior- the mean busy period of a single server in a Jackson network 
is identical to the independent server mean busy period. This is deduced by considering the limit of two 
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nearby (in time) measurements: The busy state of a server at time t and the busy state of the server at 
time t + t where r is small. If the first measurement is "busy" and the second is "idle", a busy period has 
terminated somewhere between t and t + r. In the limit of t — > this will give r times the density of busy 
period endings. On the other hand, the situation in which these to events occur is exactly what is measured 
by the correlation 

- a+Q a )g(r)a+Q a \p) = - a+Q a ) [1 + tL] a+Q a \p) + 0(r 2 ) 
However, L = L° + eW and because the first term in e vanishes, we are left with the result, 
(J|(l - a+Q a )g(T)a+Q a \p) = t(I\(1 - a+Q a )L°a+Q a \p) + 0(t 2 ), exactly the same result as for the 
unperturbed case. 

The density of the busy period endings is the inverse of the mean time between two busy periods, which 
is composed of the sum of the mean busy period and the mean idle period. Since the fraction of the of 
time that the queue is busy is also independent of the perturbation, the mean busy period must remain 
unchanged, and this result holds for all orders of the perturbation. 



6 Comparisons to Simulations 
6.1 On the fly Laplace transforms 

In order to test the accuracy of the perturbation expansion we have used simulations to compute busy-busy 
correlations for various Jackson networks. We simulate a Jackson network and monitor the queue states. 
We perform the Laplace transforms on the fly by considering the product of the instantaneous server state 
of queue a with the exponential averaged server state of queue f3. The cumulative average of these products 
are equivalent to the Laplace transform of the correlation functions at equilibrium. The simulation used is 
event driven so that it is useful to describe how these averages can be performed during an event driven 
simulation. 

If we term the "busy" indicator of the server of queue r\ at time t by b v , We consider 

/>oo 

B {u,T)= / bp(T-t)e- ut db (6.1) 
Jo 

and 



C aP {u,T)= / b a (T)b {T - t)e~ ut dt = b a (T)B p (uj,T) (6.2) 
Jo 

so that C aP {uj) =< C a0 (uj,T) > T =< b a (T)Bp(cj,T) > T . 

So that Bp(T)can be obtained by the exponential averaging defined by 

,A 

B(T + A) = e-" A [B(T ) + / e^ T b(T + r)dr] (6.3) 

Jo 

Furthermore, the busy status of the servers remain fixed between each packet arrival (or departure) event, 
thus the exponential average is very easy to evaluate in an event driven simulation. Suppose that there is 
no change in the busy state of either queue between times To and T\. Then, 



fl _ p-w(Ti-To)-, 

B^Tx) = e-^- T ^B(T ) + [ — 



bp{T +) (6.4) 

UJ 
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Averages can be obtained by integrating C a p{u),T) over time and normalizing. Consider the contribution 
to the integral of C a p(uj, T) between the times To and T\ = Tq + A: 

/ C a/3 (uj,T + T)dT = / b a (T +T)B fi (u,T +T)dT (6.5) 
Jo Jo 

but since b a , bp are fixed during this period, 

b a (T + r)Bp{u,T + r)dr = b a (T +) [ [e- ur Bp(T a +) + (1 ~ & ^ b p (T +)]dT = 

Jo w 

U) OJ 

This gives us two equations to update at each event: 



B P {T X ) = e-" A B(T )+ {1 C UA) bp(T + e) 

OJ 

1 dtC a[ ](t) = J ° dtC a0 {t) + b a (T o +)Bp{T o +) ^~^ " A) + 

b a (T Q+ )b p (T 0+ ) [A _{l-e-^) ] (g6) 

(the uj dependence is suppressed for brevity). 

Normalizing to ^ J dtC a ff(uj,t) gives an estimate of < C a p(uJ,T) >t and thus of C a p{uj). 



6.2 Simulation results 

We simulated a set of networks, for an arbitrary choice of p (defined by equation (|3. 12[) ). we measured 
off-diagonal correlations for various values of u> and subtracted the th order perturbation term p 0,9 13 from 
the measured values. The data obtained is presented in two ways: 

1. For networks that only differ by the value of the perturbation, we plot the value of the subtracted 
correlation, normalized by the perturbation SLp <a — rp^ a p,p , as a function of uj. We expect that the 
data collapse to the same curve as long as first order perturbation theory is accurate. 

2. In order to allow for data collapse of different families of networks, the subtracted values were normalized 
by the computed perturbation (equation [5101) Ppl 2 ^^ - ^pJi-^ H ^p^i'-^) + ES rj BE 'l> leaving us 
with values were plotted with respect to SLp <a — rp^ a p,p . Accurate results are reflected by a straight 
line collapse. 



A family of 2X2 matrices was tested: r^j — ^ ^ 

we scanned the range with p <G [0.05, 0.1, 0.2, 0.3, 0.32 
in figure [TJ 



, fixing {pi, P2 } = {0.3, 0.7} and {m,^} = {0.3, 0.2}, 
and plotted the results ( scaled by the perturbation) 
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Data collapse for various values of p 
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Dependence on perturbation for various Omega's 
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Figure 1: Type (1) data collapse for the Laplace transform of the (subtracted) correlation function, scaled 
by the perturbation value, the lines marked by 2p represent < 6261 > correlations, and the lines marked 
by p represent < 6162 > correlations. The deviations from the first order result (marked "theory") are due 
to higher order contributions. Error bars mark one RMS of statistical deviation ( based on splitting the 
simulation run into sub runs of 10,000 seconds). The total simulation is approximately 2,560,000 seconds , 
or 1M events on queue 1 and about 1.7M events on queue 2. The inset shows type (2) data collapse- dividing 
out the (unsealed) perturbation result and plotting the result versus the perturbation amplitude. A 45° 
straight line (marked "Theory") indicates exact matching to first order theory. Error bars are suppressed for 
clarity. 



While a good fit is observed for all values of the perturbation and of the frequency, systematic deviations, 
probably due to higher order corrections are observable. Fits of similar quality were obtained for larger 
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networks (data not shown). 



7 Conclusions 

Having described the utility of the second quantization point of view it might be considered natural to 
point out other queueing quantities that could be calculated by various manipulations of the technique; For 
example, different correlation functions such as the queue depth correlator and other moments of the queue 
depth can be obtained. Similarly, obtaining higher order terms in the perturbation expansion appears to be 
a technically non trivial extension. Indeed, a diagrammatic expansion for managing higher order terms may 
be formulatedjTTj. The dynamics of other queueing models besides M/M/l also seem to be approachable in 
this technique (as pointed out above, the M/M/oo queue and network can be written using QM harmonic- 
oscillator ladder operators). However, instead of dwelling on various extensions of the formalism (which will 
be explored elsewhere), the author would like to point out some shortcomings of the technique, in other 
words- what is lacking in the second quantization formalism. 

The reason a "reaction diffusion" approach is useful for queues is due to the fact the the questions asked 
in this paper do not relate directly to the "FIFO" nature of the problem. The particles were viewed as 
indistinguishable, and the entire information about a queue state was encoded in the number of particles 
waiting at the queue. The information about the order of the particles waiting in the queue is lost in this 
representation. Thus, questions that relate to the experience of a specific particle (e.g. waiting times) cannot 
be formulated. Issues such as "jitter" ( how the inter-arrival time of two particles at the destination is related 
to the "inter-injection" time at the source) cannot be addressed either. 

This is somewhat analogous to the difference between "Eulerian" and "Lagrangian" views in fluid 
dynamics [12]. In the Lagrangian approach one considers the trajectories of tagged fluid parcels that are 
advected with the flow. The Eulerian approach examines the behavior at fixed positions in space and con- 
siders the density and velocity of the fluid flowing through these positions. It could be argued that while the 
Eulerian approach is technically more understood, the Lagrangian approach captures various aspects that 
are very difficult to introduce in the Eulerian point of view (e.g. the invariance under Galilean invariance 
that removes the "sweeping" effect that masks various correlation functions). 

It is difficult to ask "Lagrangian" type questions it the operator formalism. Massey attempted an extension 
in |6j, but in is unclear how to use the results obtained there. In the framework of the asymmetric simple 
exclusion process (ASEP) it is possible to examine the behavior of "tracer" particles [13], but the application 
to queues is not obvious. 

A Single queue green's function 

The propagator is defined as the operator g(t) — e tc , and the Green's function is its Laplace transform 
g(oj) = We are interested in a "normal ordered form" for the Green's function in which all the Q 

operators are on the right and the a + are on the left ( so that g |0) is a simple calculation). Following [7], we 
start with 



1 



1 



(A.l) 



lo — fj,(l — a + )(Q — z) 



uj - fi{Q — p + a+p - 1) + p(a+Q - 1) 




The authors of [7] note that go(uj) = u-/j,(Q-p+a+ p-i) ls easn y normal ordered by the ansatz: 
!}xa+) (i-xQ/p) (f° r a t° be determined value of x) , which can be seen by directly applying 
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[u - (i{Q - p + a+p - 1)] {1 J xa+) {1 _* Q/p) = 

= [ u -^x-p-l)-f] {1 _l a+) {1 J Q/p) +f[l- xQ/p] JT ^ 7M 



[u-rtx-p-1)- ^] 1 1 + ^. (A.2) 

a; (1 — xa + ) (1 — xQ/p) x 



If we demand that x(uj) solves 



x 

then the above inversion implies that we can write g~o(tv) as 



[u - p(x - p -l)-E£} = 0, (A.3) 



A , . x/{pp) 1 . 

^- (l-xa^d-xQ/py (A - 4) 



Given <fo( w ) the authors of [7] perturbatively iterate for G: 

1 



g 1 + /x(a+Q - 1) 

5o - gon(a + Q ~ 1)50 + gop(a + Q - l)5oM(a + Q - l)3o 



noting that 



(a+Q - l)5o(« + Q - 1) = 10X01 ^^ (1 _^Q /p) l°)(°l = C 1 - « + QW(^) 

we retrieve 

<?M = So - ^.9o(« + Q - l)ffo = go - ^ (T^(a + - 1) (ir^y = 



x/(PP) 1 a-- / (MP) (n + D — '\\ 3- 

(l-aro+) (1-xQ/p) p-x (l-aro+) ^ ^ > (1-xC 



Q/p) 

. (1 _ _^_ (a+Q _ 1)} 1 (A.5) 



(1 — xa+) p — ir " (1 — xQ/p) 

We are left with solving equation I A.3I for a;(w): 
_ (Wm-p+i)±VWah-p+i) 2 -4|0 

• T ± — 2 

in order to decide the sign of the root, we note that for convergence we must have \x\ < 1. We further 
note that for to = 

x± = = (p + i)±\p-M = {pj 1} {p < 1} 

and since we expect convergence at this point, we must choose the negative root 



*(*,) = (^ + ^ + 1 w(^ + ^ + 1 ) 2 - 4 ^ (A . 6) 
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B Single queue busy-busy correlation 

We want to evaluate (I\Qg(uj)a + \p). The idea is to push all the Q's to the right hand side and all the a + 's 
to the left. 

<TOMa+|,) = ±<I|0__ (1 - -£-(« + Q - i))_i__ +| p ) (B.l) 

- f P {i\[Q + n^Ki - A( a+ « - i))^ + 

= /fp^ll^ + (l-i)]([(l-x) ! (p-a;) + a+ (p-J)(l-2:))l^) 

= IT^yK(l-a ; j C (p- a; )) + [1+ (p-x)(l-x))\P) 

p(l-x)(l+x)-x(l-x) lf\p\ 



Hp{l-x){p-x) 



[izp(\ — x)(p — x) 

An alternative form: 



(p + px-x). (B.2) 



/itp(l — ar) W 

can be derived by setting g — 1 _ 1 xa + g and then noting that Qg = Qg + xg. 
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